function handles = pareto2cube(P, side, alpha, linecolor, pointcolor)
%PARETO2CUBE Draw 3D pareto surface.
%   PARETO2CUBE(P) constructs a 3D pareto surface from a set
%   of pareto points P. 
%   PARETO2CUBE(P, SIDE), where SIDE is 'inside' or 'outside'
%   controls from which direction the cube is seen.
%   PARETO2CUBE(P, SIDE, ALPHA, LINECOLOR, POINTCOLOR) additionally
%   specifies the face alpha, and point- and line colors.
%   If omitted, ALPHA = 1, LINECOLOR = 'None' and
%   POINTCOLOR = 'None'

if nargin < 2
  side = 'outside';
else
  if ~strcmp(side, 'outside') & ~strcmp(side, 'inside')
    error('Side must be ''inside'' or ''outside''');
  end
end

if nargin < 3
  alpha = 1;
end

if nargin < 4
  linecolor = 'None';
end

if nargin < 5
  pointcolor = 'None';
end

handles = [];
texts = [];

% Find unique values, record forward and back mappings, and
% extend lookup tables for scaling
for i = 1:3
  [L{i}, I{i}, J{i}] = unique(P(:, i));
  EL{i} = [L{i}; max(L{i}); max(L{i})];
  PA{i} = [];
  EA{i} = [];
end

% Initialize unscaled cube
A = zeros(length(L{1}), length(L{2}), length(L{3}));

% Construct unscaled pareto-dominated accumulation cube
for p = 1:size(P, 1)
  A(J{1}(p):end, J{2}(p):end, J{3}(p):end) = 1;
end 

caxis([0 1])
colormap gray
hold on

% Create patches, for each dimension separately
for i = [1 2 3]
  % Find dimensions orthogonal to projection dimension
  D1 = min(rem(i,3), rem(i+1,3))+1;
  D2 = max(rem(i,3), rem(i+1,3))+1;
  D3 = i;

  % Construct meshes
  [M1, M2] = meshgrid(1:length(L{D1}), 1:length(L{D2}));
  M3 = max(max(sum(A, i))) - squeeze(sum(A, i))' + 1;

  % Get linear index
  E = 1:numel(M3);

  % Construct patches
  PA{D1} = [EL{D1}(M1(E)), EL{D1}(M1(E))  , EL{D1}(M1(E)+1), EL{D1}(M1(E)+1)];
  PA{D2} = [EL{D2}(M2(E)), EL{D2}(M2(E)+1), EL{D2}(M2(E)+1), EL{D2}(M2(E))  ];
  PA{D3} = [EL{D3}(M3(E)), EL{D3}(M3(E))  , EL{D3}(M3(E))  , EL{D3}(M3(E))  ]; 

  % Draw patches
  if alpha ~= 0
    h = patch(PA{1}', PA{2}', PA{3}', 0);
    set(h, 'FaceColor', repmat(i*0.25, [1 3]));
    if alpha ~= 1
      set(h, 'FaceAlpha', alpha);
      set(h, 'EdgeColor', 'None');
    else
      set(h, 'EdgeColor', repmat(i*0.25, [1 3]));
    end

    handles = [handles h];
  end

  % Construct edges
  IDX = find(M3 ~= [M3(2:end, :); M3(1, :)]);

  EA{D1} = [EA{D1}; EL{D1}(M1(IDX)), EL{D1}(M1(IDX)+1)];
  EA{D2} = [EA{D2}; EL{D2}(M2(IDX)+1), EL{D2}(M2(IDX)+1)];
  EA{D3} = [EA{D3}; EL{D3}(M3(IDX)), EL{D3}(M3(IDX))];

  IDX = IDX((rem(IDX-1, size(M3, 1))+1)<size(M3, 1));
  EA{D1} = [EA{D1}; EL{D1}(M1(IDX)), EL{D1}(M1(IDX)+1)];
  EA{D2} = [EA{D2}; EL{D2}(M2(IDX)+1), EL{D2}(M2(IDX)+1)];
  EA{D3} = [EA{D3}; EL{D3}(M3(IDX+1)), EL{D3}(M3(IDX+1))];

  IDX = find(M3 ~= [M3(:, 2:end), M3(:, 1)]);
  EA{D1} = [EA{D1}; EL{D1}(M1(IDX)+1), EL{D1}(M1(IDX)+1)];
  EA{D2} = [EA{D2}; EL{D2}(M2(IDX)), EL{D2}(M2(IDX)+1)];
  EA{D3} = [EA{D3}; EL{D3}(M3(IDX)), EL{D3}(M3(IDX))];

  IDX = IDX(IDX<=(numel(M3)-size(M3, 1)));
  EA{D1} = [EA{D1}; EL{D1}(M1(IDX)+1), EL{D1}(M1(IDX)+1)];
  EA{D2} = [EA{D2}; EL{D2}(M2(IDX)), EL{D2}(M2(IDX)+1)];
  EA{D3} = [EA{D3}; EL{D3}(M3(IDX+size(M3, 1))), EL{D3}(M3(IDX+size(M3, 1)))];
end

% Draw edges
if ~isstr(linecolor) | ~strcmpi(linecolor, 'None')
  h = line(EA{1}', EA{2}', EA{3}', 'Color', linecolor);
  handles = [handles h'];
end

% Draw pareto points
if ~isstr(pointcolor) | ~strcmpi(pointcolor, 'None')
  h = scatter3(P(:, 1), P(:, 2), P(:, 3), 40, pointcolor, 'filled');
  handles = [handles h];
end

hold off

% View mode
if strcmp(side,'inside')
  view(135, 30);
else
  set(gca, 'ZDir', 'reverse');
  view(-45, 30)
end

% Create axes
set(gca, 'XGrid', 'on');
set(gca, 'YGrid', 'on');
set(gca, 'ZGrid', 'on');

h = title('Pareto-dominated volume');
texts = [texts h];

h = xlabel('Processing time (s)');
texts = [texts h];
h = ylabel('Area (mm^2)');
texts = [texts h];
h = zlabel('Energy (J)');
texts = [texts h];

handles = [texts handles];

set(texts, 'FontSize', 18);
